#### setting environment ####
require(Matrix)
require(emIRT)
require(RColorBrewer)
color <- brewer.pal(8, "Pastel2")

account.label <- c("Asahi", "Asahi (social news)", "Asahi (sports)", "Chugoku", 
                   "Chunichi", "Chunichi (social news)", "Hokkaido (culture)", "Hokkaido (sports)", 
                   "Hokkaido", "Kahoku", "Mainichi (political news)", "Mainichi", "Mainichi (social news)", 
                   "Nikkei", "Nikkei (culture)", "Nikkei (sports)", "Nikkei (social news)", 
                   "Nishinippon", "Nishinippon (social news)", "Sankei", "Tokyo (political news)", 
                   "Yomiuri (culture)", "Yomiuri (social news)", "Yomiuri (political news)", 
                   "Yomiuri (sports)", "Yomiuri")

#### estimates by the unidimensional factor analysis model ####
# estimation results of the factor analysis models
load("FA_result.Rdata")

## newspapers' ideal points
# posterior draws
Study2.ideal.points.draws <- t(rbind(FA.result$mcmc[[1]][, 1:10], 
                                     FA.result$mcmc[[2]][, 1:10], 
                                     FA.result$mcmc[[3]][, 1:10])) * 10
# point estimates (posterior means)
Study2.ideal.points.mean <- rowMeans(Study2.ideal.points.draws)

#### Twitter network data ####
## list of Diet members' accounts
MP.list <- read.csv("MP_list.csv")

## follower-followee relationship between politically active Twitter users and Diet members
load("follower_followee_matrix.Rdata")

## follower IDs of the newspapers' accounts
load("newspaper_follower_anonymous_id.Rdata")

#### politically active Twitter users' ideal points ####
n.followers <- nrow(follower.followee.matrix)
n.followers

n.followees <- ncol(follower.followee.matrix)
n.followees

followee.name <- colnames(follower.followee.matrix)

## initial values
# initial values for Diet members' ideal points (LDP = 1, JCP = -1, others = 0)
w.inits <- matrix(NA, ncol(follower.followee.matrix), 1)
for (i in 1:ncol(follower.followee.matrix)) {
  w.inits[i, 1] <- ifelse(followee.name[i] %in% MP.list$name[MP.list$party == "LDP"], 1, 
                          ifelse(followee.name[i] %in% MP.list$name[MP.list$party == "JCP"], -1, 0))
}
emIRT.starts <- list(alpha = matrix(0, n.followees, 1), beta = matrix(0, n.followers, 1), 
                     w = w.inits, theta = matrix(0, n.followers, 1), gamma = as.matrix(1, 1, 1))

## hyperparameters of prior distributions
emIRT.priors <- list(alpha = list(mu = matrix(0, 1, 1), sigma = matrix(100, 1, 1)), 
                     beta = list(mu = matrix(0, 1, 1), sigma = matrix(100, 1, 1)), 
                     w = list(mu = matrix(0, 1, 1), sigma = matrix(100, 1, 1)), 
                     theta = list(mu = matrix(0, 1, 1), sigma = matrix(1, 1, 1)), 
                     gamma = list(mu = matrix(1, 1, 1), sigma = matrix(1, 1, 1)))

## estimate the parameters of the network ideal point model
# CAUTION: it may take a long time to complete this procedure depending on the computational environment
# you may skip this procedure and load a .Rdata file containing the estimation results
emIRT.fit <- networkIRT(.y = as.matrix(follower.followee.matrix) * 1, 
                        .starts = emIRT.starts, .priors = emIRT.priors, 
                        .control = list(verbose = TRUE, threads = 4))
# save(emIRT.fit, file = "emIRT_fit.Rdata")
# load("emIRT_fit.Rdata")

#### average of politically active followers' ideal points for each newspaper account ####
## compute average and density
newspaper.follower.ideal.point.mean <- newspaper.follower.ideal.point.SE <- rep(NA, length(newspaper.follower.anonymous.id))
newspaper.follower.ideal.point.density <- list()
for (i in 1:length(newspaper.follower.anonymous.id)) {
  newspaper.follower.ideal.point.mean[i] <- mean(emIRT.fit$means$theta[newspaper.follower.anonymous.id[[i]]])
  newspaper.follower.ideal.point.SE[i] <- sd(emIRT.fit$means$theta[newspaper.follower.anonymous.id[[i]]]) / 
    sqrt(length(emIRT.fit$means$theta[newspaper.follower.anonymous.id[[i]]]))
  newspaper.follower.ideal.point.density[[i]] <- density(emIRT.fit$means$theta[newspaper.follower.anonymous.id[[i]]])
}
names(newspaper.follower.ideal.point.mean) <- names(newspaper.follower.ideal.point.SE) <- account.label
newspaper.follower.ideal.point.order <- order(newspaper.follower.ideal.point.mean)

# average
round(newspaper.follower.ideal.point.mean, 3)

#### illustrate the results ####
# newspapers' id corresponding the estimation results of newspapers' ideal points
newspaper.id <- c(1, 1, 1, 2, 3, 3, 4, 4, 4, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 9, 3, 10, 10, 10, 10, 10)
# indicator for local newspaper accounts
local <- c(0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0)
# indicator for accounts handling "soft" topics
hard <- c(1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1)
# order of plot (plot "soft" accounts prior to "hard" accounts)
plot.order <- order(hard)

## draw a scatter plot
cairo_pdf("Figure_3_b.pdf", 
          width = 3, height = 3, pointsize = 7, family = "Helvetica")
par(mar = c(4, 4, 2.5, 2.5))
plot(NULL, NULL, type = "n", 
     xlim = c(-1, 1.5), ylim = c(-0.5, 0.4), 
     main = "(b) Following Twitter accounts", 
     xlab = "Newspapers' Ideal Points", ylab = "Average of Followers' Ideal Points")
points(Study2.ideal.points.mean[newspaper.id][plot.order], 
       newspaper.follower.ideal.point.mean[plot.order], 
       pch = ifelse(local[plot.order] == 1, 23, 21), col = "gray40", 
       bg = ifelse(hard[plot.order] == 0, "white", "gray40"))
dev.off()

## draw the distribution of politically active followers' ideal points
cairo_pdf("Figure_4.pdf", 
          width = 6, height = 7.5, pointsize = 10, family = "Helvetica")
layout(matrix(1:length(newspaper.follower.anonymous.id), length(newspaper.follower.anonymous.id), 1))
par(mar = c(0, 0, 0, 0), oma = c(4, 0, 0, 0), lwd = 0.5)
for (i in 1:length(newspaper.follower.anonymous.id)) {
  j <- newspaper.follower.ideal.point.order[i]
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(-2.3, 1.5), ylim = c(0, 3), 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  abline(v = seq(-1.5, 1.5, 0.5), lty = 3, col = "gray80")
  polygon(newspaper.follower.ideal.point.density[[j]]$x, 
          newspaper.follower.ideal.point.density[[j]]$y, border = NA, col = color[3])
  points(newspaper.follower.ideal.point.mean[j], 0.5, pch = 25, col = NA, bg = "black")
  segments(newspaper.follower.ideal.point.mean[j] + newspaper.follower.ideal.point.SE[j] * qnorm(0.025), 0, 
           newspaper.follower.ideal.point.mean[j] + newspaper.follower.ideal.point.SE[j] * qnorm(0.975), 0)
  text(-2.4, 1.5, paste0(account.label[j], "\n", 
                         length(newspaper.follower.anonymous.id[[j]]), " / ", 
                         newspaper.n.followers[j]), pos = 4)
}
axis(1, at = seq(-1.5, 1.5, 0.5), outer = TRUE, line = 0.5, lwd = 0.5)
mtext("Followers' Ideal Points", side = 1, at = 0, line = 3, cex = 0.8)
dev.off()